Relativistic electron flux decay and recovery

relative role of EMIC waves, whistler-mode waves, and plasmasheet injections

Authors

Affiliation

Zijin Zhang

University of California, Los Angeles

Anton Artemyev

University of California, Los Angeles

Motivation

We observed three strong electron precipitation which lasted for 3 hours from ELFIN missions.

We wanted to investigate the cause of the precipitation.

Import all the libraries

# To make quarto render holoview properly, somehow you can not use `output: false` in the cell options
import logging

import altair as alt
import holoviews as hv
import hvplot.pandas  # noqa
import hvplot.xarray  # noqa
import numpy as np
import pandas as pd
import panel as pn
import pdpipe as pdp
import proplot as pplt
import pyspedas
import xarray as xr
from icecream import ic
from plasmapy.formulary import gyrofrequency
from pyspedas import tinterpol
from pytplot import get_data, options, split_vec, store_data, tplot, tplot_options

from utilities import *

logger = logging.getLogger()

pplt.rc["svg.fonttype"] = "none"

Paramaters of this notebook

PLOT = True

Overview of the April 17 2021 event

Get satellite location data

Note: Location data from SSCWebServices is of low resolution (1 min). Resolution can be improved by using the data from the satellite server directly.

sats = ["elfina", "elfinb", "arase", "mms1", "goes16", "goes17", "noaa15", "noaa19"]
trange = ["2021-04-17T00:00:00", "2021-04-17T12:00:00"]

Update spacepy toolbox (updating may take a while)

import spacepy.toolbox

spacepy.toolbox.update(QDomni=True)

Retrieve data from SSCWebServices and calculate L-Shell and MLT

import spacepy
from spacepy import coordinates, irbempy, time, toolbox
from sscws.sscws import SscWs


def process_data(data, backend="IRBEM", extMag="T01STORM"):
    IRBEM_RE = 6371.2

    sat = data["Id"]
    time = pd.to_datetime(data["Time"]).tz_convert(None)
    x = data["Coordinates"][0]["X"] / IRBEM_RE
    y = data["Coordinates"][0]["Y"] / IRBEM_RE
    z = data["Coordinates"][0]["Z"] / IRBEM_RE
    coords = np.array([x, y, z]).T

    t = spacepy.time.Ticktock(data["Time"], "ISO")
    loci = spacepy.coordinates.Coords(coords, "GSE", carsph="car", use_irbem=True)
    Lm = spacepy.irbempy.get_Lm(t, loci, 90, extMag=extMag)
    # magequator=spacepy.irbempy.find_magequator(t,loci)

    return pd.DataFrame(
        {
            "Id": sat,
            "Time": pd.Series(time),
            "X": pd.Series(x),
            "Y": pd.Series(y),
            "Z": pd.Series(z),
            "Lshell": pd.Series(Lm["Lm"].flatten()).abs(),
            "MltE": pd.Series(Lm["MLT"].flatten()),
        },
    )


ssc = SscWs()
result = ssc.get_locations(sats, trange)

# Use a list comprehension to process data and store results in a list
data_frames = {data["Id"]: process_data(data, extMag="T89") for data in result["Data"]}

# Concatenate the data frames into a single DataFrame
df = pd.concat(data_frames.values())

pline = pdp.PdPipeline(
    [
        pdp.df["Angle"] << pdp.df["MltE"] * 15 / 360 * 2 * np.pi,
        pdp.df["EquatorX"]
        << pdp.df["Lshell"] * (pdp.df["MltE"] * 15 / 360 * 2 * np.pi).map(np.cos),
        pdp.df["EquatorY"]
        << pdp.df["Lshell"] * (pdp.df["MltE"] * 15 / 360 * 2 * np.pi).map(np.sin),
    ]
)

df = pline.apply(df)
for _sat, _df in data_frames.items():
    data_frames.update({_sat: pline.SetIndex("Time").apply(_df)})

Save dataframe and dataframes for later use

import pickle

df.to_pickle("data/orbit_overview.pkl")

with open("data/orbit_data_frames.pickle", "wb") as handle:
    pickle.dump(data_frames, handle, protocol=pickle.HIGHEST_PROTOCOL)

Read dataframe and dataframes from files

with open("data/orbit_data_frames.pickle", "rb") as handle:
    data_frames = pickle.load(handle)

df = pd.read_pickle("data/orbit_overview.pkl")

Inspection of the satellite orbit and sketch plot

Overview orbit plot of all satellites using hvplot

# df.hvplot("Time", "Lshell", by="Id", subplots=True) * df.hvplot("Time", "MltE", by="Id", subplots=True)

# plot individual satellites
xlim = (-15, 15)
ylim = (-15, 15)
sat_plot = df.hvplot.scatter(
    x="EquatorX",
    y="EquatorY",
    by="Id",
    subplots=True,
    aspect=1,
    xlim=xlim,
    ylim=ylim,
    hover_cols=["Time"],
)
earth = hv.Ellipse(0, 0, 2)
ellipses = hv.Overlay(
    [hv.Ellipse(0, 0, s).opts(color="gray", alpha=0.2) for s in range(6, 38, 4)]
)

# sat_plot.opts( )
sat_plot * earth * ellipses

Interactive overview plot

# adjust mms1 time to be the same as elfinb
df.loc[df["Id"] == "mms1", "Time"] = df.loc[df["Id"] == "elfinb", "Time"]

time_plot = df.hvplot.scatter(
    x="EquatorX", y="EquatorY", by="Id", groupby="Time", aspect=1
)
time_plot.opts(aspect=1, xlim=xlim, ylim=ylim)

dt_range = pn.widgets.DateRangeSlider(
    name="Date Range Slider",
    start=df.Time.unique()[0],
    end=df.Time.unique()[-1],
    format="%H:%M:%S",
)

dfi = df.interactive()
timerange_plot = dfi[
    (dfi.Time >= dt_range.param.value_start) & (dfi.Time <= dt_range.param.value_end)
].hvplot.line(x="EquatorX", y="EquatorY", by="Id", aspect=1, xlim=xlim, ylim=ylim)
timerange_plot.opts()
pn.Row(time_plot * earth * ellipses, timerange_plot * earth * ellipses)

Overview orbit plot of all satellites using proplot

groups = df.groupby("Id")
data = pd.concat({name: group.drop(columns="Id") for name, group in groups}, axis=1)

fig = pplt.figure(refwidth=5)
ax = fig.subplot(proj="polar")

ax.plot(
    data.loc[:, (sats, "MltE")] * 15 / 360 * 2 * np.pi,
    data.loc[:, (sats, "Lshell")],
    labels=sats,
)

ax.format(rlim=(0, 15))
ax.legend(loc="r", label="Orbits", ncols=1, frame=False)

Plot of satellite orbits with highlight for interested time intervals (publication-quality)

Code

fig = pplt.figure(refwidth=5)
ax = fig.subplot(proj="polar")

color_name = "Qual1"
color_index = -1
alpha = 0.618


def ax_point(time, text=True, scatter=True, marker="*"):
    time = pd.to_datetime(time)
    iloc_idx = data_frames[sat].index.get_indexer([time], method="nearest")[0]
    text_data = data_frames[sat].iloc[iloc_idx]
    rotation = text_data.Angle * 180 / np.pi
    if scatter:
        ax.scatter(
            text_data.Angle,
            text_data.Lshell,
            marker=marker,
            facecolors="None",
            edgecolors=(color_name, color_index),
        )
    if text:
        ax.text(
            text_data.Angle,
            text_data.Lshell + 1,
            time.strftime("%H:%M"),
            fontsize="small",
            rotation=rotation,
            rotation_mode="anchor",
            horizontalalignment="center",
            verticalalignment="center",
        )

Plot ELFIN-B orbit

sat = "elfinb"
tranges = [
    ["2021-04-17 02:42:00", "2021-04-17 02:46:00"],
    ["2021-04-17 04:14:00", "2021-04-17 04:18:00"],
    ["2021-04-17 05:47:00", "2021-04-17 05:51:00"],
]
color_index += 1


# create a new DataFrame containing all the time ranges
# reindex the DataFrame with the new index to fill the missing rows with NaN values
# plot the reindexed DataFrame
df = pd.concat([data_frames[sat].loc[slice(*trange)] for trange in tranges]).reindex(
    data_frames[sat].index
)

ax.plot(
    "Angle",
    "Lshell",
    data=df,
    label="ELFIN-B",
    color=(color_name, color_index),
    alpha=alpha,
)

for time in np.array(tranges)[:, 0]:
    ax_point(time)

for time in np.array(tranges)[:, 1]:
    ax_point(time, text=False, marker="D")

# plot highlighted time ranges of interest
tranges = [
    ["2021-04-17 02:43:00", "2021-04-17 02:44:00"],
    ["2021-04-17 04:15:00", "2021-04-17 04:15:00"],
    ["2021-04-17 05:48:00", "2021-04-17 05:49:00"],
]

df = pd.concat([data_frames[sat].loc[slice(*trange)] for trange in tranges]).reindex(
    data_frames[sat].index
)

ax.plot(
    "Angle",
    "Lshell",
    data=df,
    label="",
    color=(color_name, color_index),
    lw=4,
)

Plot ELFIN-A orbit (Optional)

""" This is commented out because of the low time resolution of the data from SSCWeb.
"""

if False:
    sat = "elfina"
    tranges = [
        ["2021-04-17 05:15:00", "2021-04-17 05:20:00"],
    ]
    color_index += 1

    # create a new DataFrame containing all the time ranges
    # reindex the DataFrame with the new index to fill the missing rows with NaN values
    # plot the reindexed DataFrame
    df = pd.concat(
        [data_frames[sat].loc[slice(*trange)] for trange in tranges]
    ).reindex(data_frames[sat].index)

    ax.plot(
        "Angle",
        "Lshell",
        data=df,
        label="ELFIN-B",
        color=(color_name, color_index),
        alpha=alpha,
    )

    for time in np.array(tranges)[:, 0]:
        ax_point(time)

    for time in np.array(tranges)[:, 1]:
        ax_point(time, text=False, marker="D")

    # plot highlighted time ranges of interest
    tranges = [
        ["2021-04-17 05:16:00", "2021-04-17 05:17:00"],
    ]

    df = pd.concat(
        [data_frames[sat].loc[slice(*trange)] for trange in tranges]
    ).reindex(data_frames[sat].index)

    ax.plot(
        "Angle",
        "Lshell",
        data=df,
        label="",
        color=(color_name, color_index),
        lw=4,
    )

Plot NOAA-19 orbit

sat = "noaa19"
color_index += 1
tranges = [
    ["2021-04-17 01:47:00", "2021-04-17 01:52:00"],
    ["2021-04-17 03:30:00", "2021-04-17 03:36:00"],
]

df = pd.concat([data_frames[sat].loc[slice(*trange)] for trange in tranges]).reindex(
    data_frames[sat].index
)

ax.plot(
    "Angle",
    "Lshell",
    data=df,
    label="NOAA-19",
    color=(color_name, color_index),
    alpha=alpha,
)

for time in np.array(tranges)[:, 0]:
    ax_point(time)

for time in np.array(tranges)[:, 1]:
    ax_point(time, text=False, marker="D")

tranges = [
    ["2021-04-17 01:49:00", "2021-04-17 01:51:00"],
    ["2021-04-17 03:33:00", "2021-04-17 03:35:00"],
]

df = pd.concat([data_frames[sat].loc[slice(*trange)] for trange in tranges]).reindex(
    data_frames[sat].index
)

ax.plot(
    "Angle",
    "Lshell",
    data=df,
    label="",
    color=(color_name, color_index),
    lw=4,
)

Plot NOAA-15 orbit

sat = "noaa15"
color_index += 1

tranges = [
    ["2021-04-17 01:15:00", "2021-04-17 01:20:00"],
    ["2021-04-17 02:58:00", "2021-04-17 03:03:00"],
    # ["2021-04-17 04:39:00", "2021-04-17 04:44:00"],
]

df = pd.concat([data_frames[sat].loc[slice(*trange)] for trange in tranges]).reindex(
    data_frames[sat].index
)

ax.plot(
    "Angle",
    "Lshell",
    data=df,
    label="NOAA-15",
    color=(color_name, color_index),
    alpha=alpha,
)

for time in np.array(tranges)[:, 0]:
    ax_point(time)

for time in np.array(tranges)[:, 1]:
    ax_point(time, text=False, marker="D")

tranges = [
    ["2021-04-17 01:17:00", "2021-04-17 01:19:00"],
    ["2021-04-17 02:59:00", "2021-04-17 03:02:00"],
    # ["2021-04-17 04:41:00", "2021-04-17 04:43:00"],
]

df = pd.concat([data_frames[sat].loc[slice(*trange)] for trange in tranges]).reindex(
    data_frames[sat].index
)

ax.plot(
    "Angle",
    "Lshell",
    data=df,
    label="",
    color=(color_name, color_index),
    lw=4,
)

Plot GOES-16 orbit

sat = "goes16"
color_index += 1

tranges = [
    ["2021-04-17 00:00:00", "2021-04-17 12:00:00"],
]
df = pd.concat([data_frames[sat].loc[slice(*trange)] for trange in tranges]).reindex(
    data_frames[sat].index
)
ax.plot(
    "Angle",
    "Lshell",
    data=df,
    label="GOES-16",
    color=(color_name, color_index),
    alpha=alpha,
)
for time in np.array(tranges)[:, 0]:
    ax_point(time)

for time in np.array(tranges)[:, 1]:
    ax_point(time, text=False, marker="D")

tranges = [
    ["2021-04-17 02:00:00", "2021-04-17 02:30:00"],  # ion injection
    ["2021-04-17 04:00:00", "2021-04-17 04:30:00"],  # electron injection
    ["2021-04-17 07:00:00", "2021-04-17 07:30:00"],  # electron injection
]

df = pd.concat([data_frames[sat].loc[slice(*trange)] for trange in tranges]).reindex(
    data_frames[sat].index
)

ax.plot(
    "Angle",
    "Lshell",
    data=df,
    label="",
    color=(color_name, color_index),
    lw=4,
)

for time in np.array(tranges).flatten():
    ax_point(time, scatter=False)

Plot MMS-1 orbit

# ! NOTE: mms1 obtained Lm does not agree with the one provided (doesn't matter for this plot)
sat = "mms1"
color_index += 1


tranges = [
    ["2021-04-17 00:00:00", "2021-04-17 12:00:00"],
]
df = pd.concat([data_frames[sat].loc[slice(*trange)] for trange in tranges]).reindex(
    data_frames[sat].index
)

ax.plot(
    "Angle",
    "Lshell",
    data=df,
    label="MMS-1",
    color=(color_name, color_index),
    alpha=alpha,
)

for time in np.array(tranges)[:, 0]:
    ax_point(time)

for time in np.array(tranges)[:, 1]:
    ax_point(time, text=False, marker="D")


tranges = [
    ["2021-04-17 04:30:00", "2021-04-17 04:45:00"],  # electron decrease
    ["2021-04-17 08:00:00", "2021-04-17 08:15:00"],  # electron injection
]

df = pd.concat([data_frames[sat].loc[slice(*trange)] for trange in tranges]).reindex(
    data_frames[sat].index
)

ax.plot(
    "Angle",
    "Lshell",
    data=df,
    label="",
    color=(color_name, color_index),
    lw=4,
)

for time in np.array(tranges).flatten():
    ax_point(time, scatter=False)

Plot ARASE orbit

sat = "arase"
color_index += 1


tranges = [
    ["2021-04-17 00:00:00", "2021-04-17 12:00:00"],
]
df = pd.concat([data_frames[sat].loc[slice(*trange)] for trange in tranges]).reindex(
    data_frames[sat].index
)

ax.plot(
    "Angle",
    "Lshell",
    data=df,
    label="ARASE",
    color=(color_name, color_index),
    alpha=alpha,
)
for time in np.array(tranges)[:, 0]:
    ax_point(time)

for time in np.array(tranges)[:, 1]:
    ax_point(time, text=False, marker="D")

tranges = [
    ["2021-04-17 01:00:00", "2021-04-17 01:15:00"],  # electron injection
    ["2021-04-17 07:10:00", "2021-04-17 07:25:00"],  # electron injection
]

df = pd.concat([data_frames[sat].loc[slice(*trange)] for trange in tranges]).reindex(
    data_frames[sat].index
)

ax.plot(
    "Angle",
    "Lshell",
    data=df,
    label="",
    color=(color_name, color_index),
    lw=4,
)

for time in np.array(tranges).flatten():
    ax_point(time, scatter=False)

Plot the Earth

# create an array of angles from 0 to 2*pi
r = np.ones(1000)
theta = np.linspace(0, 2 * np.pi, 1000)

# plot the circle
ax.plot(theta, r, color="black", lw=0.5)

# fill the left half of the circle with black
ax.fill_between(theta, r, where=theta < np.pi / 2, color="black")
ax.fill_between(theta, r, where=theta > np.pi * 3 / 2, color="black")

Format the plot

def degree_to_hour(degrees, pos):
    """A custom formatting function to convert degrees to hours"""
    hour = (degrees / 2 / np.pi) * 24
    return f"{int(hour)}"


ax.format(
    title="Multi-mission orbits 2021-04-17 00:00 to 12:00 UTC",
    thetaformatter="func",
    thetaformatter_kw={"func": degree_to_hour},
    rlim=(0, 15),
    rlabelpos=135,
)

ax.legend(loc="r", ncols=1, frame=False, title="Orbits")

fig

An overview of the mission orbits recorded on April 17, 2021, from 00:00 to 12:00 UTC. The orbits of distinct missions are projected onto the MLT and L-Shell plane, designated with different colors; star markers denote the orbit initiation, squares indicate their termination, and time annotations are provided near the periods of interest. ELFIN-B’s trajectory is demonstrated in three time intervals: 02:42-02:46, 04:14-04:18, and 05:47-05:51. NOAA-19’s trajectory is plotted from 01:47-01:52 and 03:30-03:36, while NOAA-15’s is displayed from 01:15-01:20 and 02:58-03:03. The trajectories of GOES, MMS, and ARASE span the entirety of the 12-hour interval from 00:00 to 12:00.

Detailed analysis of the April 17 2021 event

ELFIN observations of EMIC driven precipitations during >3hour interval

Two ELFIN satellites observation of EMIC-driven precipitation (where the precipitation flux surpasses the trapped flux in high-energy channels) over an interval exceeding three hours, from 02:42 to 05:53. The locations are projected in proximity to the L-Shell and MLT, using the Tsyganenko (1989) magnetic field model. Panels (a), (b), and (d) encapsulate data from ELFIN-B, while panel (c) features observations from ELFIN-A.

ARASE Part

ARASE observation of whistler wave

Process data from the ERG satellite: retrieve data, interpolate to the same time as the spectra, and export tplot variables

Retrieve data from the ARASE satellite

trange = ["2021-04-17 00:00:00", "2021-04-17 12:00:00"]
pyspedas.erg.pwe_ofa(trange=trange, time_clip=True, ror=False)
pyspedas.erg.mgf(trange=trange, time_clip=True, ror=False)
pyspedas.erg.orb(trange=trange, time_clip=True, ror=False)
pyspedas.erg.orb(
    trange=trange, level="l3", model="t89", time_clip=True, ror=False
)  # Get equatorial magnetic field data

# Interpolate the spacecraft position and magnetic field data to the same time as the spectra
split_vec(
    "erg_orb_l2_pos_rmlatmlt", new_name="erg_orb_l2_pos_", suffix=["r", "mlat", "mlt"]
)

tinterpol("erg_mgf_l2_magt_8sec", "erg_pwe_ofa_l2_spec_B_spectra_132")
tinterpol("erg_orb_l2_pos_Lm", "erg_pwe_ofa_l2_spec_B_spectra_132")
tinterpol("erg_orb_l2_pos_mlat", "erg_pwe_ofa_l2_spec_B_spectra_132")
tinterpol("erg_orb_l3_pos_beq_t89", "erg_pwe_ofa_l2_spec_B_spectra_132")

# Get the data from pytplot variables
erg_mgf_magt_itrp = get_data("erg_mgf_l2_magt_8sec-itrp")
pwe_spec = get_data("erg_pwe_ofa_l2_spec_B_spectra_132")
pwe_spec_xr = get_data("erg_pwe_ofa_l2_spec_B_spectra_132", xarray=True)
orb_Lm_itrp = get_data("erg_orb_l2_pos_Lm-itrp", xarray=True)
orb_mlat_itrp = get_data("erg_orb_l2_pos_mlat-itrp", xarray=True)
orb_beq_itrp = get_data("erg_orb_l3_pos_beq_t89-itrp", xarray=True)

Overview of ARASE observations

# Plot the PWE data with the gyrofrequency data and magnetic field data

# Store gyrofrequency data with the PWE data
mag_t = erg_mgf_magt_itrp.y * u.nT

omega_e = gyrofrequency(mag_t, "e-", to_hz=True).to("kHz")
upper_bound = 0.1 * gyrofrequency(mag_t, "e-", to_hz=True).to("kHz")
lower_bound = 0.5 * gyrofrequency(mag_t, "e-", to_hz=True).to("kHz")

store_data("omega_e", data={"x": erg_mgf_magt_itrp.times, "y": omega_e})
store_data("omega_e01", data={"x": erg_mgf_magt_itrp.times, "y": lower_bound})
store_data("omega_e05", data={"x": erg_mgf_magt_itrp.times, "y": upper_bound})

# convert omega_e to a DataArray
omega_e_da = xr.DataArray(omega_e, coords=[pwe_spec_xr.time], dims=["time"])

# store the combined data
store_data(
    "combined_data",
    data=["omega_e", "omega_e01", "omega_e05", "erg_pwe_ofa_l2_spec_B_spectra_132"],
)

options("combined_data", "spec", True)
options("omega_e", "yrange", [pwe_spec.v[0], pwe_spec.v[-1]])

tplot(["combined_data"])
tplot(["erg_orb_l2_pos_Lm", "erg_orb_l3_pos_blocal_t89", "erg_orb_l3_pos_beq_t89"])

# tplot(["combined_data"],save_png='../figures/erg_pwe_spec')

22-Jun-23 13:06:01: combined_data does not contain coordinates for spectrogram plotting.  Continuing...

Get the whistler waves normalized mean frequencies from Arase mission for calculating diffusion coefficients.

Code

# filter the data
lower_bound = 0.1
upper_bound = 0.5
step = 0.025

interested_freq = lambda x: (x.spec_bins < upper_bound * omega_e_da) & (
    x.spec_bins > lower_bound * omega_e_da
)
pwe_spec_filtered = pwe_spec_xr.where(interested_freq, 1e-5)

# only keep the time when the maximum intensity is at the interested frequency (whistler mode)
interested_time = pwe_spec_xr.argmax(dim="v_dim") == pwe_spec_filtered.argmax(
    dim="v_dim"
)
pwe_spec_temp = pwe_spec_xr.where(interested_time)

pwe_spec_filtered = pwe_spec_temp.where(interested_freq)


def summed_intensity(_lower_bound, _upper_bound):
    interested_freq = lambda x: (x.spec_bins < _upper_bound * omega_e_da) & (
        x.spec_bins > _lower_bound * omega_e_da
    )
    # TODO: better way to acculmulate this intensity for uneven frequency bins
    pwe_spec_freq = (
        pwe_spec_temp.where(interested_freq).mean(dim=["v_dim"])
        * (_upper_bound - _lower_bound)
        * omega_e_da
        * 1000  # convert to Hz
    )

    rel_freq = xr.DataArray(
        [(_upper_bound + _lower_bound) / 2],
        dims=["relative_frequency"],
        attrs={"long_name": "Relative Frequency", "units": r"$\omega_{ce}$"},
    )
    pwe_spec_freq.attrs["long_name"] = "Wave intensity"
    pwe_spec_freq.attrs["units"] = r"$pT^2$"

    return pwe_spec_freq.expand_dims(
        dim={"relative_frequency": rel_freq}
    ).assign_coords({"relative_frequency": rel_freq})


lower_bounds = np.arange(lower_bound, upper_bound, step)
upper_bounds = np.arange(lower_bound + step, upper_bound + step, step)
da = xr.concat(
    [
        summed_intensity(_lower_bound, _upper_bound)
        for _lower_bound, _upper_bound in zip(lower_bounds, upper_bounds)
    ],
    dim="relative_frequency",
)

avg_da = (
    da.assign_coords({"Lm": ("time", orb_Lm_itrp.data[:, 0])})
    .where(lambda x: (x.Lm < 6.5) & (x.Lm > 4.5))
    .mean(dim="time")
)
relative_frequency = avg_da.idxmax().values
print(f"Normalized mean frequencies of whistler waves: {relative_frequency}")

# Plot the summed intensity of the PWE data between 0.1 and 0.5 omega_e
# fig, axs = pplt.subplots()
# axs.plot(avg_da)

Normalized mean frequencies of whistler waves: 0.31249999999999994

Whistler mode waves observed by ARASE

# Calculate omega_e_max
omega_e_max = relative_frequency * gyrofrequency(mag_t, "e-", to_hz=True).to("kHz")

# store the filtered data
store_data(
    "pwe_spec_filtered",
    data={"x": pwe_spec.times, "y": pwe_spec_filtered, "v": pwe_spec.v},
)

# Store omega_e_max data
store_data("omega_e_max", data={"x": erg_mgf_magt_itrp.times, "y": omega_e_max})

# Store the combined data
combined_data = ["omega_e_max", "pwe_spec_filtered"]
store_data("combined_data", data=combined_data)

# Set options
options("combined_data", "spec", True)
options("pwe_spec_filtered", "spec", True)
options("pwe_spec_filtered", "ylog", True)
options("pwe_spec_filtered", "yrange", [pwe_spec.v[0], pwe_spec.v[-1]])
options("pwe_spec_filtered", "ytitle", "Frequency")
options("pwe_spec_filtered", "ysubtitle", "[kHz]")
options("pwe_spec_filtered", "zlog", True)
options("pwe_spec_filtered", "zrange", [1e-4, 1e2])
options("pwe_spec_filtered", "ztitle", "Intensity")
options("pwe_spec_filtered", "zsubtitle", "[pT^2/Hz]")

# Plot the combined_data
tplot("combined_data")

21-Jun-23 21:17:22: combined_data does not contain coordinates for spectrogram plotting.  Continuing...

ARASE electron and proton observations

Retrieve ARASE data

from pyspedas import erg

trange = ["2021-04-17T00", "2021-04-17T12"]

vars = {
    "erg_mep_ele": "erg_mepe_l2_omniflux_FEDO",
    "erg_mep_pro": "erg_mepi_l2_omniflux_FODO",
    "erg_mgf": "erg_mgf_l2_mag_8sec_gsm",
    "erg_mlt": "erg_orb_l2_pos_mlt",
    "erg_l": "erg_orb_l2_pos_Lm",
}

erg.mepe(trange=trange, time_clip=True, ror=False)
erg.mepi_nml(trange=trange, time_clip=True, ror=False)
erg.mgf(trange=trange, time_clip=True, ror=False, varformat="*mag*")
erg.orb(trange=trange, time_clip=True, ror=False, varformat="*pos*")

split_vec(
    "erg_orb_l2_pos_rmlatmlt", new_name="erg_orb_l2_pos_", suffix=["r", "mlat", "mlt"]
)

split_vec("erg_orb_l2_pos_Lm", suffix=["_90", "_60", "_30"])

Overview of ARASE observations

tplot(list(vars.values()))

Get ARASE data and update the metadata

ds = {key: get_data(value, xarray=True) for (key, value) in vars.items()}


for key in ds.keys():
    match key:
        case "erg_mep_ele":
            ds[key].attrs.update(
                {
                    "long_name": "ARASE electron",
                    "units": r"$1/s/cm^2/sr/keV$",
                    "labels": [f"{erg:0.0f} keV" for erg in ds[key].spec_bins.values],
                }
            )
            # Replace all occurrences of 0 with NaN
            ds[key] = ds[key].where(ds[key] != 0)
        case "erg_mep_pro":
            ds[key].attrs.update(
                {
                    "long_name": "ARASE proton",
                    "units": r"$1/s/cm^2/sr/keV$",
                    "labels": [f"{erg:0.0f} keV" for erg in ds[key].spec_bins.values],
                }
            )
            # Replace all occurrences of 0 with NaN
            ds[key] = ds[key].where(ds[key] != 0)
        case "erg_mgf":
            ds[key].attrs.update(
                {
                    "long_name": "ARASE B",
                    "units": "nT",
                    "labels": [r"$B_x$", r"$B_y$", r"$B_z$"],
                }
            )
        case "erg_mlt":
            ds[key].attrs.update({"long_name": "MLT"})
        case "erg_l":
            ds[key].attrs.update({"long_name": "L shell"})
            # Replace all occurrences of 0 with NaN
            ds[key] = ds[key].where(ds[key] != 0)

Plot ARASE electron and proton observations

def plot_erg_mep(type=None, save=False):
    fname = "arase_mep"

    if type == "poster":
        n = 3
        fname = f"{fname}_poster"
    else:
        n = 4

    fig, axs = pplt.subplots(
        ncols=1, nrows=n, refaspect=n, refwidth=2.5 * n, sharey=False
    )

    axs[0].plot(
        ds["erg_mep_ele"][:, ::2],
        cycle="darkred",
        labels=ds["erg_mep_ele"].attrs["labels"][::2],
    )
    axs[1].plot(
        ds["erg_mep_pro"][:, 3:-3:2],
        cycle="darkblue",
        labels=ds["erg_mep_pro"].attrs["labels"][3:-3:2],
    )

    axs[n - 1].plot(ds["erg_mlt"], color="black")
    axs[n - 1].alty(color="red", label=ds["erg_l"].attrs["long_name"]).plot(
        ds["erg_l"], color="red"
    )

    axs[0].format(yscale="log", yformatter="log", ymin=1e3)
    axs[0].legend(loc="right", ncols=1, frame=False)
    axs[1].format(yscale="log", yformatter="log")
    axs[1].legend(loc="right", ncols=1, frame=False)

    if not type == "poster":
        axs[2].plot(
            ds["erg_mgf"],
            cycle=["black", "red", "blue"],
            labels=ds["erg_mgf"].attrs["labels"],
        )
        axs[2].format(ylim=(-3e2, 3e2))
        axs[2].legend(loc="right", ncols=1, frame=False)

    if save:
        fig.savefig(f"../figures/{fname}.svg")
        fig.savefig(f"../figures/{fname}.pdf")

    return fig


fig = plot_erg_mep(type="poster", save=True)

ARASE electron and proton observations (~10 keV to ~100 keV). Strong electron injection are visible at the beginning of EMIC-driven electron precipitations and at the end of interval.

Investigate changing fo high energy fluxes over time from ARASE

Load the data

trange = ["2021-04-16", "2021-04-18"]

pyspedas.erg.xep(trange=trange, time_clip=True, ror=False)
pyspedas.erg.orb(trange=trange, time_clip=True, ror=False)

# xep: 'erg_xep_l2_FEDO_SSD'
# orb: ['erg_orb_l2_pos_llr',
#  'erg_orb_l2_pos_gse',
#  'erg_orb_l2_pos_gsm',
#  'erg_orb_l2_pos_sm',
#  'erg_orb_l2_pos_rmlatmlt',
#  'erg_orb_l2_pos_eq',
#  'erg_orb_l2_pos_iono_north',
#  'erg_orb_l2_pos_iono_south',
#  'erg_orb_l2_pos_blocal',
#  'erg_orb_l2_pos_blocal_mag',
#  'erg_orb_l2_pos_beq',
#  'erg_orb_l2_pos_beq_mag',
#  'erg_orb_l2_pos_Lm',
#  'erg_orb_l2_vel_gse',
#  'erg_orb_l2_vel_gsm',
#  'erg_orb_l2_vel_sm',
#  'erg_orb_l2_spn_num',
#  'erg_orb_l2_man_prep_flag',
#  'erg_orb_l2_man_on_flag',
#  'erg_orb_l2_eclipse_flag']

Plot the overview of the high energy fluxes for ARASE

tplot("erg_xep_l2_FEDO_SSD")
tplot("erg_orb_l2_pos_Lm")

Get data from tplot variables

split_vec(
    "erg_orb_l2_pos_rmlatmlt", new_name="erg_orb_l2_pos_", suffix=["r", "mlat", "mlt"]
)

erg_orb_Lm = get_data("erg_orb_l2_pos_Lm", xarray=True)
erg_orb_mlat = get_data("erg_orb_l2_pos_mlat", xarray=True)
erg_orb_mlt = get_data("erg_orb_l2_pos_mlt", xarray=True)
erg_xep = get_data("erg_xep_l2_FEDO_SSD", xarray=True)

Filter and group the data

def time_group(time):
    """Group data points by time"""

    # Set a threshold time interval for grouping data
    time_threshold = pd.Timedelta(minutes=30)

    # Compute the time intervals between each data point
    time_delta = time.diff()

    # Create a boolean mask for identifying the start of each group
    group_start_mask = (time_delta > time_threshold) | time_delta.isnull()

    # Assign a group ID to each data point based on the start of each group
    return group_start_mask.cumsum()


erg_orb_pipe = pdp.PdPipeline(
    [
        # only keep data where Lm > 5 and 10 > mlat > -10
        pdp.keep_rows_where["erg_orb_l2_pos_Lm"] > 4.5,
        pdp.keep_rows_where["erg_orb_l2_pos_Lm"] < 5.5,
        pdp.keep_rows_where["erg_orb_l2_pos_mlat"] > -10,
        pdp.keep_rows_where["erg_orb_l2_pos_mlat"] < 10,
        # group data by time
        pdp.ColByFrameFunc("group_id", lambda df: time_group(df.index.to_series())),
    ]
)

erg_orb_df = xr.merge([erg_orb_Lm[:, 0], erg_orb_mlat, erg_orb_mlt]).to_dataframe()
erg_orb_df_filtered = erg_orb_pipe.apply(erg_orb_df)

Plot the overview of the ERG orbit during the period of interest

erg_orb_df.hvplot.scatter(
    x="erg_orb_l2_pos_mlt", y="erg_orb_l2_pos_mlat"
) * erg_orb_df_filtered.hvplot.scatter(x="erg_orb_l2_pos_mlt", y="erg_orb_l2_pos_mlat")
erg_orb_df.hvplot(
    x="time", y=["erg_orb_l2_pos_Lm", "erg_orb_l2_pos_mlat", "erg_orb_l2_pos_mlt"]
) * erg_orb_df_filtered.hvplot.scatter(
    x="time", y=["erg_orb_l2_pos_Lm", "erg_orb_l2_pos_mlat", "erg_orb_l2_pos_mlt"]
)

Group erg_xep data by orbit groups

erg_orb_tgroup = (
    erg_orb_df_filtered.reset_index()
    .groupby("group_id")
    .agg({"time": ["min", "max"]})
    .to_numpy()
)

erg_xep_df = erg_xep.to_dataframe()

for tr_min, tr_max in erg_orb_tgroup:
    erg_xep_df.loc[tr_min:tr_max, "trange"] = "{}, {}".format(
        tr_min.astype("datetime64[h]"), tr_max.astype("datetime64[h]")
    )

Code

vatt = erg_xep.attrs["CDF"]["VATT"]

chart01 = (
    alt.Chart(erg_xep_df.dropna())
    .mark_line(point=True)
    .encode(
        x=alt.X("spec_bins", title="Energy [keV]"),
        y=alt.Y(
            "mean(erg_xep_l2_FEDO_SSD)",
            scale=alt.Scale(type="log"),
            title=f"{vatt['LABLAXIS']} [{vatt['UNITS']}]",
        ),
        color="trange",
    )
    .properties(title=f"{vatt['CATDESC']}")
)

chart02 = (
    alt.Chart(erg_xep_df.dropna())
    .mark_line(point=True)
    .encode(
        x=alt.X("spec_bins", title="Energy [keV]"),
        y=alt.Y(
            "mean(erg_xep_l2_FEDO_SSD)", title=f"{vatt['LABLAXIS']} [{vatt['UNITS']}]"
        ),
        color="trange",
    )
    .properties(title=f"{vatt['CATDESC']}")
)

display(chart01, chart02)

MMS electron and proton observations

Retrieve MMS data

from pyspedas import mms, time_clip

probe = "1"
trange = ["2021-04-17T00", "2021-04-17T12"]

vars = {
    "mms1_feeps_ele": "mms1_epd_feeps_srvy_l2_electron_intensity_omni_spin",
    "mms1_feeps_proton": "mms1_epd_eis_srvy_l2_extof_proton_flux_omni_spin",
    "mms1_fgm_bvec": "mms1_fgm_b_gsm_srvy_l2_bvec",
    "mms1_mlt": "mms1_mec_mlt",
    "mms1_l_dipole": "mms1_mec_l_dipole",
}

mms.feeps(
    trange=trange, probe=probe, varformat="*electron*"
)  # BUG:`varnames` doesn't work
mms.eis(
    trange=trange, probe=probe
)  # BUG:`varnames` doesn't work ; `varformat` doesn't load spin data
mms.fgm(trange=trange, probe=probe, varformat="*gsm*")  # BUG:`varnames` doesn't work
mms.mec(trange=trange, probe=probe, varformat="*mec*")  # BUG:`varnames` doesn't work

for var in vars.values():
    time_clip(var, trange[0], trange[1], suffix="")

Overview of MMS observations

tplot(list(vars.values()))

22-Jun-23 11:40:19: /Users/zijin/mambaforge/envs/cool_space_science/lib/python3.10/site-packages/pytplot/MPLPlotter/specplot.py:102: RuntimeWarning: divide by zero encountered in log10
  zdata = np.log10(out_values)

Get MMS data into xarray and update the metadata

ds = {key: get_data(value, xarray=True) for (key, value) in vars.items()}

for key in ds.keys():
    match key.split("_", 1)[1]:
        case "feeps_ele" | "feeps_proton":
            ds[key].attrs.update(
                {
                    "long_name": ds[key].attrs["plot_options"]["yaxis_opt"][
                        "axis_label"
                    ],
                    "units": ds[key].attrs["plot_options"]["zaxis_opt"]["axis_label"],
                    "labels": [f"{erg:0.0f} keV" for erg in ds[key].spec_bins.values],
                }
            )
        case "fgm_bvec":
            ds[key].attrs.update(
                {
                    "long_name": "MMS1 B",
                    "units": "nT",
                    "labels": [r"$B_x$", r"$B_y$", r"$B_z$"],
                }
            )
        case "mlt":
            ds[key].attrs.update({"long_name": "MLT"})
        case "l_dipole":
            ds[key].attrs.update({"long_name": "L shell"})

Plot MMS electron and proton observations

def plot_mms_epd(type=None, save=False):
    fname = "mms_epd"
    if type == "poster":
        n = 3
        fname = f"{fname}_poster"
    else:
        n = 4

    fig, axs = pplt.subplots(
        ncols=1, nrows=n, refaspect=n, refwidth=2.5 * n, sharey=False
    )
    axs[0].plot(
        ds["mms1_feeps_ele"][:, ::2],
        cycle="darkred_r",
        labels=ds["mms1_feeps_ele"].attrs["labels"][::2],
    )
    axs[1].plot(
        ds["mms1_feeps_proton"][:, :-1:],
        cycle="darkblue_r",
        labels=ds["mms1_feeps_proton"].attrs["labels"][:-1],
    )

    axs[n - 1].plot(ds["mms1_mlt"], color="black")
    axs[n - 1].alty(color="red", label=ds["mms1_l_dipole"].attrs["long_name"]).plot(
        ds["mms1_l_dipole"], color="red"
    )

    axs[0].format(yscale="log", yformatter="log", ymin=1e1)
    axs[0].legend(loc="right", ncols=1, frame=False)
    axs[1].format(yscale="log", yformatter="log", ymin=1e0)
    axs[1].legend(loc="right", ncols=1, frame=False)

    if not type == "poster":
        axs[2].plot(
            ds["mms1_fgm_bvec"],
            cycle=["black", "red", "blue"],
            labels=ds["mms1_fgm_bvec"].attrs["labels"],
        )
        axs[2].format(ylim=(-3e2, 3e2))
        axs[2].legend(loc="right", ncols=1, frame=False)

    if save:
        fig.savefig(f"../figures/{fname}.svg")
        fig.savefig(f"../figures/{fname}.pdf")

    return fig


# plot_mms_epd(type='poster', save=True)
fig = plot_mms_epd()

MMS electron and proton observations (~50 keV to ~500 keV). Localized decrease of electron fluxes is notable around ELFIN observations of EMIC driven precipitations and strong electron injection signatures are visible at the end of EMIC-driven electron precipitations.

POES observations of proton precipitation

Overview of the POES from THEMIS website

POES observations of proton precipitation

logger.disabled = True

probe = "noaa15"
time_slices = [
    ["2021-04-17 01:15:00", "2021-04-17 01:20:00"],
    ["2021-04-17 02:57:00", "2021-04-17 03:03:00"],
    ["2021-04-17 04:38:00", "2021-04-17 04:44:00"],
]

[
    plot_poes_mep_pro(probe=probe, trange=time_slice, save=True)
    for time_slice in time_slices
]

probe = "noaa19"
time_slices = [
    ["2021-04-17 01:47:00", "2021-04-17 01:52:00"],
    ["2021-04-17 03:30:00", "2021-04-17 03:36:00"],
    ["2021-04-17 05:12:00", "2021-04-17 05:18:00"],
]

[
    plot_poes_mep_pro(probe=probe, trange=time_slice, save=True)
    for time_slice in time_slices
]

logger.disabled = False

POES observations of proton precipitation confirms the L-shell localization of EMIC generation region between L~4.5 and 6. At L>6 POES shows strong curvature scattering of protons indicating on isotropic equatorial distribution and no EMIC wave generation”

GOES observations of electron and proton

Overview of GOES from THEMIS website

Particle data is from Space Environment In Situ Suite (SEISS) instrument.

Retrieve GOES particle data

trange = ["2021-04-17T00:00:00", "2021-04-17T12:00:00"]

g16_pflux, g16_eflux = goes_part_prod(trange, probe="16")
g17_pflux, g17_eflux = goes_part_prod(trange, probe="17")

g16_mag = geos_mag_prod(trange, probe="16")
g17_mag = geos_mag_prod(trange, probe="17")
# pyspedas.goes.orbit(trange=trange, probe=probe, time_clip=True, prefix=f'g{probe}_')

Plot GOES electron and proton observations

def plot_goes(type=None):
    if type == "poster":
        n = 4
    else:
        n = 6

    fig, axs = pplt.subplots(
        ncols=1, nrows=n, refaspect=n, refwidth=2.5 * n, sharey=False
    )

    axs[0].plot(
        g16_eflux.sum(dim="telescopes"),
        cycle="darkred_r",
    )
    axs[1].plot(
        g16_pflux.sum(dim="telescopes"),
        cycle="darkblue_r",
    )
    axs[2].plot(
        g17_eflux.sum(dim="telescopes"),
        cycle="darkred_r",
    )
    axs[3].plot(
        g17_pflux.sum(dim="telescopes"),
        cycle="darkblue_r",
    )

    axs[0].format(
        yscale="log",
        yformatter="log",
        ylim=(1e-1, 1e6),
    )
    axs[1].format(
        yscale="log",
        yformatter="log",
        ylim=(1e-2, 1e5),
    )
    axs[2].format(
        yscale="log",
        yformatter="log",
        ylim=(1e-1, 1e6),
    )
    axs[3].format(
        yscale="log",
        yformatter="log",
        ylim=(1e-2, 1e5),
    )

    axs[0].legend(loc="right", ncols=1, frame=False)
    axs[1].legend(loc="right", ncols=1, frame=False)
    axs[2].legend(loc="right", ncols=1, frame=False)
    axs[3].legend(loc="right", ncols=1, frame=False)

    if not type == "poster":
        axs[4].plot(g16_mag)
        axs[5].plot(g17_mag)
        axs[4].legend(loc="right", ncols=1, frame=False)
        axs[5].legend(loc="right", ncols=1, frame=False)

    if type == "poster":
        fig.savefig(f"../figures/goes_mpsh_poster.svg")
        fig.savefig(f"../figures/goes_mpsh_poster.pdf")
    elif type == "publication":
        fig.savefig(f"../figures/goes_mpsh.svg")
        fig.savefig(f"../figures/goes_mpsh.pdf")
    return fig


# fig = plot_goes(type='poster', save=True)
fig = plot_goes()

GOES-R electron and proton observations (70 keV to ~1 MeV) from two geostationary operatioinal satellites. Ion injections are visible before the time when ELFIN observed strong electron precipitation. Series of strong electron injections observed around noon after drifting from the midnight.